1. Field of the Invention
The present invention relates to an image processing apparatus for performing an image processing, for example, including a dynamic range compression processing for adding a high frequency component of an original image to the original image (dynamic range change processing), an image processing system, an image processing method, and a memory medium for storing processing steps to perform the processing in a computer readable manner.
2. Related Background Art
In recent years, with an advancement of digital technique, for example, a method is executed which comprises: digitizing a photographed image obtained by X-ray photographing (hereinafter also referred to as “X-ray image”), performing an image processing to the digitized image; and displaying the image on a monitor apparatus, or outputting the image onto an X-ray diagnosis film.
Here, the X-ray image is constituted of an image area in which it is easy to transmit X-rays, and an image area in which it is very difficult to transmit the X-rays. For example, the X-ray image of a chest (lungs) is constituted of a lung image area in which it is easy to transmit the X-rays, and a mediastinum image area in which it is very difficult to transmit the X-rays. In this manner, since the dynamic range in which pixel values are present (hereinafter also referred to as “density range”) is much broadened in the X-ray image, it has been difficult to obtain an X-ray image in which both the lung and mediastinum image areas far different in X-ray transmittance from each other can simultaneously be observed.
Therefore, in order to avoid the above-described problem, the following methods 1 to 5 of a dynamic range compression processing (hereinafter also referred to as “DRC processing”) have been proposed.
Method 1:
In method 1 described in “SPIE Vol. 626 Medicine XIV/PACSIV (1986)” or the like, assuming that the pixel value of the processed image is “SD”, the pixel value of the original image (input original image) is “Sorg”, and the pixel value of the low frequency component of the original image (pixel value of a smoothed image) is “SUS”, equation (1) is represented with constants A, B, C (e.g., A=3, B=0.7).SD=A[Sorg−SUS]+B[SUS]+C  (1)
Moreover, in the method 1, by changing weighting (constants A and B) of the high frequency component (first term) and low frequency component (second term) in the equation (1) and, for example, assuming A=3, B=0.7, the high frequency component is emphasized, and the dynamic range of the entire image can be compressed. This is evaluated by many radiotherapists, et al. such that the image subjected to the present processing is effective for diagnosis as compared with the image not subjected to the present processing.
Method 2:
In a method 2 described in Japanese Patent Publication No. 6-046409, assuming that the pixel value of the processed image is “SD”, the pixel value of the original image is “Sorg”, and the pixel value of the low frequency component of the original image (pixel value of the smoothed image) is “SUS”, equation (2) is represented with a monotonous decrease function f(X).SD=Sorg+f(SUS)  (2)
Also in this method 2, similarly to the above-described method 1, the dynamic range of the entire image can be compressed based on the low frequency component of the original image.
Method 3:
In a method 3 described in Japanese Patent No. 2509503, assuming that the pixel value of the processed image is “SD”, and the pixel value of the original image is “Sorg”, equation (3) is represented with the average profile Py of the Y-directional profile of the original image, and the average profile Px of the X-directional profile.SD=Sorg+F[G(Px,Py)]  (3)
Here, the property of the function F(x) in the equation (3) will be described. First, in “x>Dth”, F(0) is “0”. Moreover, in “0≦x≦Dth”, F(x) monotonously decreases with an intercept “E”, and inclination “E/Dth”, and is represented by equation (4).F[x]=E−(E/Dth)X  (4)
Moreover, the average profiles Py and Px in the equation (3) are represented by equations (5) and (6) with profiles pyi and pxi (i=1 to n).Py=(ΣPyi)/n  (5)Px=(ΣPxi)/n  (6)
Moreover, “G(Px,Py)” in the equation (3) is represented by equation (7).G(Px,Py)=max(px,py)  (7)
Therefore, in the method 3, the dynamic range is compressed when the pixel value of the low frequency component is Dth or less.
Method 4:
A method 4 similar to the method 2 described in the Japanese Patent Publication No. 6-046409 and the method 3 described in the Japanese Patent No. 2509503 is described in “Journal of Japanese Society of Radiological Technology, Vol. 45, No. 8, August, 1989, p. 1030, Anan et al”.
Assuming that the pixel value of the processed image is a “SD”, the pixel value of the original image is “Sorg”, and the average pixel value (pixel value of the smoothed image) obtained by taking the moving average of the original image with a mask size of M×M pixels is “SUS”, the method 4 is represented by equations (8) and (9) with the monotonous decrease function f(X).SD=Sorg+f(SUS)  (8)SUS=ΣSorg/M2  (9)
Moreover, the equation (8) can be changed to equation (10).                                                                         S                D                            =                            ⁢                                                (                                                            S                      org                                        -                                          S                      US                                                        )                                +                                  (                                                            f                      ⁡                                              (                                                  S                          US                                                )                                                              +                                          S                      US                                                        )                                                                                                        =                            ⁢                                                (                                                            S                      org                                        -                                          S                      US                                                        )                                +                                  f1                  ⁡                                      (                                          S                      US                                        )                                                                                                          (        10        )            
Here, the method 4 is different from the method 3 represented by the equation (3) in method of generating the low frequency component, the low frequency component is generated with one-dimensional data in the method 3, and the low frequency component is generated by two-dimensional data in the method 4.
Also in the method 4, similarly to the above-described method 3, the dynamic range is compressed when the pixel value of the low frequency component is Dth or less.
Method 5:
For a method 5 described in Japanese Patent No. 2663189, assuming that the pixel value of the processed image is a “SD”, the pixel value of the original image is “Sorg”, and the average pixel value (pixel value of the smoothed image) obtained by taking the moving average of the original image with the mask size of M×M pixels is “SUS”, equations (11) and (12) are represented with the monotonous decrease function f2(X).                                                                         S                D                            =                            ⁢                                                S                  org                                +                                  f2                  ⁡                                      (                                          S                      US                                        )                                                                                                                          =                            ⁢                                                (                                                            S                      org                                        -                                          S                      US                                                        )                                +                                  f3                  ⁡                                      (                                          S                      US                                        )                                                                                                                                          f3                ⁡                                  (                                      S                    US                                    )                                            =                            ⁢                                                f2                  ⁡                                      (                                          S                      US                                        )                                                  +                                  S                  US                                                                                        (        11        )            SUS=ΣSorg/M2  (12)
Here, the property of the function f2(x) in the equation (11) will be described. First, in “x<Dth”, f2(0) is “0”. Moreover, in “Dth≦x”, f2(x) monotonously decreases with the intercept “E”, and inclination “E/Dth”, and is represented by equation (13).f2[x]=E−(E/Dth)X  (13)
Therefore, in the method 5, the dynamic range is compressed when the low frequency component pixel value is Dth or less.
Additionally, the compression algorithm of the dynamic range in the method 5 is similar to the algorithm in the method 4 described in the “journal of Japanese Society of Radiological Technology, Vol. 45, No. 8, August, 1989, p. 1030, Anan et al”.
However, the conventional image processing method using the above-described methods 1 to 5 of the DRC processing has at least the following problems 1 and 2.
Problem 1:
For example, as shown in FIG. 14, when the original image (input original image) is subjected to the DRC processing represented by the “Journal of Japanese Society of Radiological Technology, Vol. 45, No. 8, August, 1989, p. 1030, Anan et al.”, particularly to the DRC processing including the processing (original image−smoothed image) of subtracting the pixel value SUSof the smoothed image from the pixel value Sorg of the original image as shown in the equations (10) and (11), artifacts called overshoot and undershoot are generated in the edge part of the inputted original image.
Concretely, first the cause of the overshoot and undershoot will be described with reference to FIGS. 15A, 15B and 15C.
In FIGS. 15A, 15B and 15C, solid lines show image profiles, and broken lines show coordinates. Moreover, FIG. 15A shows the profile of the edge part of the original image, FIG. 15B shows the profile of the image (smoothed image) obtained by smoothing the original image, and FIG. 15C shows the profile of the image (corresponding to the image of the high frequency component) obtained after subtracting the smoothed image shown in FIG. 15B from the original image shown in FIG. 15A.
As shown in FIG. 15C, when the original image is smoothed, the shape of the original image profile is not stored in the edge part. This causes the overshoot and undershoot.
On the other hand, in FIG. 14, the abscissa indicates an image coordinate, and the ordinate indicates an image pixel value. Moreover, in FIG. 14, (A) shows an original image profile. Here, the original image is a step-shaped image in which the pixel value increases by 450 at each 300 pixels, or an image in which the pixel value of 50 as the high frequency component is added to every 50 pixels. Moreover, in FIG. 14, (B) shows the profile of the image obtained by performing the conventional DRC processing to the original image.
Therefore, as apparent from (B) of FIG. 14, since the conventional DRC processing includes a processing of subtracting the smoothed image (SUS) of the original image from the original image (Sorg) to extract the high frequency component, for the reason described with reference to FIG. 15, the artifacts (overshoot and undershoot) are generated in the edge part of the original image.
Problem 2:
The object of the DRC processing by the above-described methods 1 and 2 is to compress the dynamic range of the pixel values of the object area in the image and obtain the image in which the entire object area can simultaneously be observed.
Therefore, for example, for the lung X-ray image, by compressing the dynamic range of the mediastinum area with a low density in the lung and mediastinum image areas constituting the lung area, it is possible to obtain the image in which both the lung and mediastinum image areas can simultaneously be observed. This is because the density of the low-density side image area is raised to a visible area by compressing the dynamic range of the low-density side image area which is not included in the visible area.
However, in the X-ray image, the low-density side image area of the object area is very low in X-ray transmittance, and usually the S/N ratio also tends to be low as compared with the (high-density side) area high in X-ray transmittance. Therefore, noises are sometimes conspicuous at a certain pixel value or less. In this case, even when the pixel value of the low-density side image area is raised to the visible area by the DRC processing, there is a problem that the effective information of the object area, particularly the information of the low frequency component such as density distribution is obstructed by the noises and cannot easily be observed.